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We investigate the importance of including quantized initial conditions in Langevin dynamics 
for adsorbates interacting with a thermal reservoir of electrons. For quadratic potentials the time 
evolution is exactly described by a classical Langevin equation and it is shown how to rigorously 
obtain quantum mechanical probabilities from the classical phase space distributions resulting from 
the dynamics. At short time scales, classical and quasiclassical initial conditions lead to wrong results 
and only correctly quantized initial conditions give a close agreement with an inherently quantum 
mechanical master equation approach. With CO on Cu(lOO) as an example, we demonstrate the 
effect for a system with ab initio frictional tensor and potential energy surfaces and show that 
quantizing the initial conditions can have a large impact on both the desorption probability and the 
distribution of molecular vibrational states. 
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I. INTRODUCTION 

Femtosecond lasers has proven a most valuable tool in 
the study of excited metallic electrons and their interac- 
tions with surface adsorbates. In Ref. [ij it was shown 
that a femtosecond laser pulse could be used to desorb 
NO from Pd(lll) and a mechanism involving multiple 
electronic excitations of the adsorbate was identified. 0^11] 
Since then, it has been demonstrated that several other 
surface reactions can be induced by femtosecond laser 
pulses and the mechanism is usually attributed to 

a direct interaction of excited (hot) metallic electrons in- 
teracting with adsorbate resonant states, although sub- 
strate heating may also contribute to reaction rates. [TT| 

A variety of theoretical models have been proposed 
to describe the interaction and resulting transfer of en- 
ergy from hot electrons to adsorbates, but a common 
conceptual picture is can be given in terms of Born- 
Oppenheimer potential energy surfaces. It is then as- 
sumed that the adsorbate propagation is governed by 
a potential energy surface Vq when the adsorbate is in 
its electronic ground state. If the adsorbate has a reso- 
nance (possibly partly occupied in the ground state), a 
hot metallic electron can transiently occupy the resonant 
state and the adsorbate dynamics will then be governed 
by a different potential energy surface Vi. Hot electrons 
can thus transfer energy to the adsorbate by inducing 
jumps between the two potential energy surfaces, i Al- 
though the lifetime of the excited electronic state on the 
adsorbate may be very short (< 1 fs), several such events 
can eventually transfer enough energy for the adsorbate 
to overcome a reaction barrier. 

The probability that a hot electron scatters inelasti- 
cally on the adsorbate and transfers a given amount of 
energy can be calculated in a local polaron mo del 1 121- 
[T5| and may be generalized to reactions resulting from 
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multiple electronic excitations. [la| However, since we are 
usually only interested in the adsorbate dynamics, it is 
often more convenient to apply open system density ma- 
trix theory. In this formalism, it is assumed that the 
femtosecond laser pulse gives rise to a hot thermalized 
distribution of electrons with a time dependent electronic 
temperature Tg. The time dependent density matrix of 
the full interacting system is then constructed and the 
electronic states are traced out resulting in a reduced 
density matrix with a diagonal that gives the probabili- 
ties that the adsorbate is in a particular state. Based on 
the Feynman- Vernon theory of influence functionals. [TtI- 
IToj it is possible to calculate the reduced density matrix 
of a Newns- Anderson type Hamiltonian in either a co- 
ordinate basis ,20j leading to Langevin dynamics or in 
a basis of vibrational eigenstates[2l| leading to a mas- 
ter equation for the vibrational eigenstates. For a har- 
monic potential with frequency ujq, the master equation 
reduces to a Fokker-Planck equation in the classical limit 
of fcsTe ^ Hu/Q and desorption probabilities can be ob- 
tained from an Arrhenius type ex pre ssion. [2^ However, 
as shown explicitly in Refs. [2Jand[2l|, the Fokker-Planck 
equation fails dramatically when the classical condition 
above is not satisfied and in general a quantum mechan- 
ical treatment of the adsorbate is needed. On the other 
hand, the coordinate representation of the reduced den- 
sity matrix results in semi-classical dynamics for the ad- 
sorbate coordinates and the quantum nature of the prob- 
lem only enters through the initial state. 

Langevin dynamics have been applied with reasonable 
success to problems involving hot electron induced sur- 
face reactions [il, HI] and to elucidate the role of non- 
adiabatic effects in general, [l^ [l^ However, the initial 
quantum state is usually neglected or treated quasiclas- 
sically. The purpose of the present work is to investigate 
the role of quantum mechanical boundary conditions and 
compare the results to those obtained with classical and 
quasiclassical initial states where only the zero point en- 
ergy is included. In particular, we will focus on the har- 
monic oscillator since, when the initial state is included 
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correctly, Langevin dynamics with a quadratic potential 
is exact to second order in perturbation theory and we 
can thus compare with a quantum mechanical master 
equation. 

The paper is organized as follows. In section |lT] we 
introduce the model Hamiltonian which constitutes the 
foundation of the calculations. The time dependent den- 
sity matrix of the harmonic oscillator is then reviewed 
and is shown to give rise to classical dynamics with quan- 
tum corrections entering only through the initial state 
which must be included by a phase space sampling pro- 
cedure. Generalizing this approach to our model Hamil- 
tonian results in Langevin dynamics with explicit expres- 
sions for the electronic friction tensor and correlations 
between fluctuating forces. In section Hill we start by an- 
alyzing the harmonic oscillator and show how to obtain 
the quantum mechanical probabilities from the classical 
phase space distribution resulting from a Langevin equa- 
tion approach. It is demonstrated that, when the ini- 
tial conditions is correctly taken into account, the results 
show excellent agreement with the master equation ap- 
proach. The comparison is then repeated for the Morse 
potential where the Langevin dynamics does not provide 
an exact description of the quantum dynamics, but which 
has the advantage of having a well defined desorption en- 
ergy. In section IIVI we consider the example of hot elec- 
tron induced desorption of CO from Cu(lOO) using ab 
initio potential energy surfaces, and perform Langevin 
dynamics with classical, quasiclassical, and quantum me- 
chanical initial conditions. In appendix]^ it is show how 
classical dynamics and the initial Wigner phase space dis- 
tribution emerges from a path integral representation of 
the time dependent reduced density matrix in a quadratic 
potential. 



II. THEORY 

A. Hamiltonian 

The Langevin dynamics with local electronic friction 
can be derived from a Newns- Anderson [2^ [2^ type 
Hamiltonian where a single adsorbate resonant state \a) 
is coupled to the adsorbate degrees of freedom Xi.[2Qii The 
resonant state is usually chosen as an eigenstate of the 
adsorbate far from the surface. Close to the surface, \a) 
becomes hybridized with metallic states and acquires a 
finite lifetime. In the electronic ground state, the reso- 
nant state has a partial (or zero) occupation and the ad- 
sorbate propagation is governed by a ground state Born- 
Oppenheimer potential energy surface Vo{xi) with a local 
minimum at x^. However, the presence of hot metallic 
electrons may give rise to a transient full occupation of 
the resonant state and the adsorbate propagation will 
then be governed by the potential energy surface ^1(2:^). 
Even though the resonant state is short lived, a transient 
occupation will perturb the system and may result in a 
transfer of energy to the adsorbate [T^. The Hamiltonian 



describing the system can then be modelled bv [l^ - [l^ 

H = H,i + Ha + Hi, (1) 

Hel = Socica + ^ Efccj' Cfe + ^ V^^clck + h.C. 
k k 

i 

Hi = {ea{xt) - eo)clca + ^ {Vak{xi) - K°fc)clcfc + h.C. 

k 

£a[x.i) = Vi[x.i) - Vo{xi) 

where c\ and c|, are creation operators for the reso- 
nant state I a) and metallic states \k) respectively and 
= Saixf), V^^. = Vak{.x^i). Conceptually, the Hamilto- 
nian describes an adsorbate with dynamics governed by 
Vo(a;i) in the electronic ground state and Vi {xi) when the 
resonant state is occupied, and the reservoir of metallic 
electrons can exchange energy with the adsorbate via the 
resonant state. The hybridization depends on the posi- 
tion of the adsorbate through Vak{xi) which become zero 
when the adsorbate is far from the surface. It should be 
noted that if Vak are constant and the ground and excited 
state potentials are quadratic with displaced minima, one 
obtains Hi = —c\ca fiXi. The coupling constants are 



then given by 



rriiUj'jxi 



where Xi is the shift in the 



minimum of the excited state potential with respect to 
the ground state minimum. 

We will impose the wide band limit in which the metal- 
lic band coupled to the adsorbate is assumed to be much 
wider than the resonance width. For a fixed position of 
the adsorbate, the density of states projected onto the 
resonance is then a Lorentzian: 



Pais) 



1 



r/2 



7r(£~£a)2 + (r/2)2' 

with the full width at half maximum given by 



r = 27r 



k 



\Vak?5{ea~ek). 



(2) 



(3) 



In these expressions both Vak and Sa and therefore pa and 
r depend parametrically on the instantaneous position of 
the adsorbate. 



B. The density matrix 

The advantage of the density matrix formalism is two- 
fold. First of all, for complicated systems one may trace 
out all irrelevant degrees of freedom from the density ma- 
trix and the resulting 'reduced' density matrix then de- 
scribes a system which can exchange energy with the en- 
vironment. Second, the density matrix formalism allows 
one to treat a statistical ensemble of states in a natural 
way. In the case of an adsorbate interacting with elec- 
trons in a metal, as described by the Hamiltonian (jlj. 
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the full density matrix can be reduced by tracing out 
the electronic degrees of freedom and the diagonal ele- 
ments of the resulting reduced density matrix then gives 
the probabilities of finding the adsorbate in a particular 
state as a function of time. 

The time dependent density matrix is 



equal to the time evolution of a classical phase space 
distribution. 
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U\ -iHt/h iHt/h 

p{t) = e ' pae ' , 



(4) 



where po is the density matrix at t = Q. As always it is 
instructive to consider a harmonic oscillator and we thus 
start by considering Hq of Eq. ([1]) with a single degree 
of freedom and a quadratic potential. In the coordinate 
representation the density matrix can then be written 



p{x,y]t) = {x\p{t)\y) ^ J dxodyo{xQ\po\yo) (5) 
X {x\e'^^°'^^\xo){yo\e'''°'^%). 

We see that the density matrix involves two propagators 
and the integrand can be viewed as a particle first be- 
ing propagated forward in time from xq to x and then 
backward in time from y to j/o- The propagator of the 
harmonic potential is well known [30| and the result for 
the diagonal elements is 

p{u;t)= J duQdpoV{uo,po) 

J./ , . . sinwi \ 
X d\u[t) — [uocosujt + Po J (6) 

where 

V(x,p) = ^ y dv{x - y/2|pok + y/2)e*f^/'' (7) 

is the Wigner distribution of an initial state described by 
the density matrix po and u = x = y. The Wigner distri- 
bution is often referred to as a quasi probability distribu- 
tion and can be interpreted as the quantum mechanical 
probability of finding a particle in the small phase space 
area dxdp. [sij This means that the expression (|6]) can be 
thought of as a sum over all initial phase space config- 
urations weighted by their probabilities and subject to 
the constraint dictated by the delta function. However, 
the constraint is equivalent to the Newtonian equations 
of motion and we can thus regard the time evolution as 
purely classical. In particular, given an initial state we 
could calculate p{u] t) by sampling all phase space and 
adding V{uq,pq) if uq and pq is classically connected to 
u{t). Furthermore, since each such classical trajectory 
will result in a well defined momentum at time t we inter- 
pret the probability of being at a given phase space point 
u{t),p{t) as being equal to V{uo,po) where (uo,po) is the 
unique point which is classically connected to {u{t),p{t)). 
The quantum nature of the particle propagating in a har- 
monic oscillator potential thus solely enters through the 
initial state specified by po- This is of course closely 
related to the well known fact that for a harmonic po- 
tential, the time evolution of the Wigner distribution is 



The Langevin equations emerge when the electronic 
degrees of freedom is traced out from the time depen- 
dent density matrix corresponding to the full Hamilto- 
nian ([T|). With a quadratic potential the result is very 
similar to ([6]) the only difference being that the coupling 
to a thermal reservoir of electrons introduces a broaden- 
ing in the delta function. Thus the time evolution can 
be thought of as classical with fluctuations that has a 
magnitude determined by the broadening. It has previ- 
ously been shown that these fluctuations can be handled 
in a statistical sensefH, [l^ and the full dynamics can 
be written in terms of classical equations of motion with 
a stochastic force (,i{t). The stochastic force is specified 
by its statistical properties which is related to the broad- 
ening of the delta function. The result is the Langevin 
equation 

where the local temperature dependent friction tensor is 
given by 

x/,(e;u)/,(e;u)^^^4^ (9) 



with 



de 



ea{u)-e dT{u) dea{u) 



du. 



(10) 



being the (dynamical) frictional force on the mode Ui. 
This result was derived in Ref. [23] for a single adsorbate 
mode and has been generalized to more than one modes 
here. It is also straightforward to extend the derivation 
to include N resonant states and the resulting friction is 
simply the sum of the N partial frictions resulting from 
each resonance. The diagonal elements of the friction ten- 
sor are strictly positive and the main contribution from 
rjij in (IS]) will be a frictional force in a direction opposite 
the velocity. In the presence of hot metallic electrons, the 
ground state potential appearing in (jS]) should actually 
be replaced by a temperature dependent renormalized 
potential Vb(u,i) Vo{ui) + F{T;Ui).^ However, the 
correction is usually so small that it can be neglected and 
we have explicitly verified this for the systems considered 
in the present work. 

In the present work we will make the Markov approx- 
imation where there is no temporal correlation of the 
fluctuating forces. The approximation is valid when the 
thermal correlation time tc ~ h/ksT is much smaller 
than the timescale of adsorbate motion, and the fluc- 
tuating force ^i(t) is a Gaussian distributed stochastic 
variable with a correlation function given by 



(6(ii)Cj(i2)) = 2r^,jkBT6iti - t2). 



(11) 
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To summarize, the Langevin equation ([8]) can be 
thought of as describing classical dynamics with stochas- 
tic fluctuations. Quantum effects enters through the ini- 
tial state of the adsorbate and can be included by run- 
ning classical trajectories with initial conditions sampled 
from a Wigner distribution of the initial state. For non- 
quadratic potentials the Langevin equation should be re- 
garded as a semiclassical approximation to the true dy- 
namics. The derivation leading to Eqs. (|8|- ((T0| is based 
on a path integral representation of the reduced density 
matrix. [20| In appendix |X] we derive Eq. ([SJ and show 
how the Wigner distribution emerges in this formalism 
and using the technique of Brandbyge et al. [20| it is 
straightforward to generalize the result to the full Hamil- 
tonian ([T]). 



C. Master equation 

If one is interested in the time dependent probability 
for the adsorbate to be be in a particular energy eigen- 
state rather than at certain position, it is more conve- 
nient to consider the reduced density matrix in a basis of 
Hamiltonian eigenstates. Taking the electronic trace of 
the Liouville equation leads to 



dpr 



dt 



■[HQ,Pr 



-Tre/ [Hi,p]. 



(12) 



where pred = Tre/(p) is the reduced density matrix and 
Tre/ is the trace over electronic states. In a basis of eigen- 
states of Hq, the diagonal elements of the reduced density 
matrix are the time dependent probabilities of finding the 
adsorbate in a particular state. The right hand side is a 
complicated functional which depends on the complete 
history of the density matrix. However, making the self 
consistent Born approximation, the Markov approxima- 
tion and neglecting the off-diagonal elements of pred leads 
to the master equation plj| 



dpn 
dt 



OO 

[pniWni-^n - P7iWn^m^ , (13) 



m=0 



where p„ = {pred)nn and Wm^n are the transition rates 
given by 



W„ 



27r 



^nF{eq){l - nF{eq'))\{q;m\Hi\q';n)\' 



9.9' 



X S{Sq ~ Eq' + Sn Srri), 



(14) 



where \q) is the eigenstates of Hf,i with eigenenergies Sq 
and npis) is the Fermi-Dirac distribution. 



III. MODEL POTENTIALS 

As shown above, zero point motion (or any other ini- 
tial quantum state) can be included in the molecular dy- 
namics by sampling all phase space and weighing each 



point according to the Wigner distribution of the initial 
state. For Langevin dynamics this can be tedious work 
since one has to run a large number of trajectories for 
each initial point in phase space to get reasonable statis- 
tics. An often used approximation to avoid phase space 
sampling is to use the classical initial conditions which 
reproduces the energy of the initial quantum state En. 
When the friction is small compared to the period of 
oscillation, one can then use a single initial phase space 
point with -Ecias (a^OiPo) = -^n- We will refer to this as the 
quasiclassical approximation. However, as will be shown 
below, this method can give rise to seriously misleading 
results for Langevin dynamics when the timescale of the 
hot electron pulse is sufficiently short. 



A. Quadratic potential 

For a quadratic potential the Langevin equation is ex- 
act within second order perturbation theory provided we 
include the initial quantum state properly. We can thus 
compare results obtained by integrating the Langevin 
equation with those obtained from a master equation 
approach (I13p and transition rates calculated from the 
Fermi golden rule expression (jl4p . In principle, the two 
approaches should be equivalent since the level of approx- 
imation is the same (Markov approximation and second 
order perturbation theory) and we can investigate the 
importance of using quasiclassical initial conditions com- 
pared to true quantum initial conditions. 

It may be surprising that the classical Langevin equa- 
tion should give the same result as the master equation 
which is inherently quantum mechanical. Furthermore, 
it may not be obvious how the probabilities p„, which 
is the basic quantity calculated within the master equa- 
tion approach, can be extracted from Langevin dynam- 
ics. However, if one has access to the Wigner distribution 
V{x,p) at a given time, it is indeed possible to calculate 
Pn since 



Pn 



= {n\p\n)^J dxdyp{x,y)(pl{x)(pniy) 
= J dudvp{u + v/2,u — v/2) 

X / dvipniu + v/2)ipn{u ~ v/2)6{v — v) 



(15) 



dudvpiu + v/2, u — v/2) 



1 



X dvdp^*n{u + i/2)ipn{u-i/2)e'P^''-^'^/'' 

= 27rh J dudpVniu,p)V{u,p), 

where Vn(u,p) is the Wigner distribution of the pure 
state density matrix /?„ — \n){n\. Integrating the 
Langevin equation gives rise to a final state classical 
phase space distribution, but since the equation of mo- 
tion for a classical phase space distribution is identical to 



5 



that of a Wigner distribution in a harmonic potential, [31 1 
we can identify the final state classical phase space dis- 
tribution with the final state Wigner distribution. 

The pure state Wig ner distributions in a quadratic po- 
tential is given by[31[ 

('— 1 "l" 

Vn{x,p) = ^-^e-^^^^P^/''°Ln{2Hix,p)/Eo), (16) 
nil 

where H{x,p) — p^/2m + ma;^a;^/2 is the classical Haniil- 
tonian, Eq = fkjj/2, and L„ is the n'th Laguerre polyno- 
mial. Since P„ is only a function of the Hamiltonian 
energy we can write 



Pn 



2nh / dEVn{E) 



dP 
dE 



(17) 



with 



= 2{-ir dEe-^/^'>K{2E/Eo)^, 
dP f 

— = J dxdpr{x,p)6{E-n{x,p)). 



(18) 



Note that the distribution dP/dE is not a true proba- 
bility distribution since it is not strictly positive, but it 
can be rigorously translated into the quantum mechan- 
ical probabilities p„. On the other hand, we can obtain 
the distribution dPn/dE associated with a particular vi- 
brational state \n) by replacing V{x,p) in Eq. (jlSp with 
Vn{x,p)- Using that dxdp = hd(pd'H/2EQ with ip being a 
phase space angle, the integral can then be evaluated to 

^ - ^-^e-^/^"K{2E/Eo). (19) 

The distributions Eq. ((T9)) are shown in Fig. [1] for the 
first four vibrational states with Eq = 0.125 eV. The 
structure of the distributions is in sharp contrast to that 
obtained in the quasiclassical (QC) approach where the 
energy is fixed at En and the energy distribution of the 
7i'th state is dPj.'^'^^/dE = S{E ~ En) with En = hw{n + 
1/2). This gives rise to completely different and and even 
negative probabilities. For example, using dP^^^^/dE = 
5{E — Eq) immediately yields po — pi = —p2 — 0.74 from 
Eq. (HZl). 

We have performed Langevin dynamics using Eqs. 
([8]) and ^ with a single mode and a linear interac- 
tion Hamiltonian: Hi = — fc\caX using the parameters 
m = 6.86 amu, hu = 0.25 eV/, Eq = 2.6 eV, F = 2.0 eV, 
and / = 8.7 eV/ A. These parameters were chosen to 
mimic the internal vibrational mode of CO adsorbed 
on Cu(lOO) considered below, but presently we will just 
think of them as a realistic set of parameters which we use 
to compare different model calculations. The adsorbate 
is initially in its ground state described by the Wigner 
distribution 

Vo{xo,Po) - -\e-^o/-l-pl/pl (20) 
Trn. 
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FIG. 1: The energy distributions given by Eq. (I19() for 
the lowest four vibrational states of a harmonic oscillator 
with zero point energy Eo — 0.125 eV. The correspond- 
ing quasiclassical distributions are deltafunctions centered at 
-Eo(2n-f 1). 



with the quantum length and momentum given by 



xq = y^h/muj, 



PQ 



(21) 



The distribution is even in both momentum and posi- 
tion and since the frictional decay is much slower than 
the vibrational time of oscillation, the final state phase 
space distribution can be assumed to be even in the initial 
phase space point. For simplicity we assume a constant 
electronic temperature at Tg = 4000 K and integrate 
the Langevin equation for t = 1 ps. For each point on 
an initial (6x6) positive phase space grid with a spac- 
ing 0.5a;Q X O.Spg, we run a large number of Langevin 
trajectories (^ 30000) and record the final state en- 
ergy. The final state energy distribution is then obtained 
by summing the distributions resulting from each initial 
phase space point dP/dE{E; xo,po) weighted by the ini- 
tial state Wigner distribution V{xo,po): 



dP{E) 
dE 



dxadpoV{xo,po) 



dP{E;xo,po) 
dE 



(22) 



In Fig. [2] we show this distribution at t = 0.1 ps and t = 
0.5 ps along with the distributions resulting from quasi- 
classical (initial phase space points with 'H(xo,po) ^ Eq) 
and classical initial condition (initial phase space point 
2^0 = Po =0). On long time scales the distributions will 
forget the initial conditions and approach a Boltzmann 
distribution at the appropriate temperature. However, 
on timescales less than a picosecond there is still plenty 
of memory of the initial state and the classical and qua- 
siclassical distributions, which start as delta functions at 
E = and E = Eq respectively, are completely wrong 
at timescales on the order of 0.1 ps. The quasiclassi- 
cal initial conditions approach the correct distribution 
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FIG. 2: The continuous energy distributions dP/dE obtained 
from Langevin dynamics with a constant Te = 4000 using 
quantum, quasiclassical, and classical boundary conditions. 
The initial quantum state is the vibrational ground state. 
Left: t = 0.1 ps. Right: t = 0.5 ps. After a while both the 
quasiclassical and classical distributions approach the quan- 
tum distribution. 



faster than the classical one since the initial state con- 
tains the right amount of energy which just needs to be 
redistributed. 

With the interaction Hamiltonian Hi = — fc\caX it is 
easy to calculate the transition rates Eq. ([T4| with the 
result: 

7r/2 f 

Wm^n ^rnS.m^n+ij-j- j ds p a{e) p a{e + fno) 
X np{e){l — np[e + huj)) 

TTp f 

+ {m + l)Sm,n-ijj^ / depa{e)pa{s - hu) 
xnFis){l-nFie-hio)). (23) 

Using the parameters above we can then integrate the 
master equation Eq. p3p and compare the probabili- 
ties pn with those obtained from the Langevin equation 
Eqs. nil) and This is shown in Fig. [3] for the 

four lowest vibrational states. As expected we see a close 
correspondence between the master equation approach 
and Langevin dynamics with correct phase space sam- 
pling. In contrast, the classical initial conditions result 
in completely wrong probabilities and the quasiclassical 
initial conditions only result in sensible probabilities after 
0.5 ps. 

It should be noted, that the quasiclassical initial con- 
ditions gives a good description of average quantities and 
the average energy (E) = ^^Pn^n is very well ap- 
proximated by the quasiclassical approach, even at short 
timescales. However, if one were to model a surface re- 
action with a barrier by a truncated harmonic potential 
the quasiclassical approach is likely to fail. For exam- 
ple, the adsorption energy of CO on Cu(lOO) is ~ 0.6 eV 
and as a simple model for hot electron induced desorp- 
tion one could use the present oscillator truncated above 
the desorption energy. This means that p2 + Pa would 
be a measure of the desorption probability and from Fig. 
[3] it is clear that for times < 0.5 ps one would severely 
miscalculate the desorption probability. 
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FIG. 3: The time dependent probabilities p„ for being in 
the vibrational state |n) obtained with the master equation 
and Langevin dynamics with three kinds of initial conditions. 
The correct quantum initial conditions are seen to give results 
nearly identical to the master equation, whereas the classical 
and quasiclassical initial conditions give wrong results. For 
small time scales the classical and quasiclassical initial con- 
ditions are not shown since the are not consistent with the 
harmonic oscillator Wigner distribution in the sense that they 
give rise to probabilities which are negative or larger than one. 



B. Morse potential 

Allthough the quadratic potential comprises a nice 
toy model for comparing Langevin dynamics with the 
master equation approach, it is not particularly well 
suited to simulate surface reactions such as desorption 
or dissociation. We will make a simple model for a 
desorption potential and modify the quadratic potential 
considered above to a one-dimensional Morse potential 
Vm{x) = I?(l-e-''^)2 with D = 0.57 eV. The parameter 
a is determined by requiring that the second derivative at 
the minimum of the well match the frequency of the har- 
monic potential considered above. A quantization of this 
potential yields five bound states with energies En and a 
continuum of free states with energies Ek = fi^k^ /2m. 

Under the influence of a thermal pulse of electrons, 
a bound state |m) can make transitions to other bound 
states \n) or to free states The transition rates can be 
calculated within second order perturbation theory and 
the result is 



W„ 



27r/2|(m|a;|n)|^ 



mn ) 

y.nF{e){l-nF{e + hhJmn)) (24) 
for bound state transitions and 



T^m^fc^ ^ i\J I /I I depaie)paie + hwmk) 



X n_F(e)(l - n_F(e -|- humk)) 



(25) 
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FIG. 4: Desorption probabilities as a function of the elec- 
tronic temperature Te calculated from the master equation 
approach and Langevin dynamics with classical quasiclassi- 
cal and quantum initial conditions. The four figures show 
the desorption probability after interaction times of 0.25 ps, 
0.5 ps, 0.75 ps, and 1.0 ps respectively. 



for transitions to free states. Here we have defined 
humi — Em — Ei. The matrix elements have been cal- 
culated previously [HI and it is now straightforward to 
integrate the master equation ([T^ . We will interpret the 
probability of being in a free state \k) at time t as the 
desorption probability. 

For a non-quadratic potential the Langevin equation is 
based on a semiclassical approximation. However, since 
the master equation l|13p is still correct within second 
order perturbation theory we can explicitly examine the 
validity of the semiclassical approximation by compar- 
ing the two approaches. Due to the lack of a classi- 
cal/quantum correspondence for the Morse potential it 
is not possible to convert the classical energy distribu- 
tion resulting from Langevin dynamics into probabilities 
of being in eigenstates of the Morse potential. Neverthe- 
less, it is natural to associate the probability of being in a 
continuum state with the probability that a classical 
trajectory results in a final state energy Ek- The initial 
quantum state is included as described above by sam- 
pling phase space and integrate weighting by the Wigner 
distribution. The Wigner distribution of the Morse po- 
tential ground state is well known [33], but since it is not 
even in the position coordinate we need to sample twice 
the phase space compared with the harmonic oscillator. 

The desorption probabilities calculated with the mas- 
ter equation and Langevin dynamics is shown in Fig. |4l 
For t — 0.25 ps, the probabilities show significant devi- 
ation signalling a breakdown of classical time evolution 
at small time scales which is expected. It is a bit more 
surprising, that the high temperature limit deviates from 
the quantum probabilities even at t = \ ps. This could 
be due a breakdown of perturbation theory at such high 
temperatures, since the effective perturbation of the sys- 




1.15 1.20 1.25 
C-0 distance [Al 



FIG. 5: Potential energy surfaces for the ground and excited 
state of CO adsorbed at a Cu(lOO) top site. The contours are 
at 0.05 eV intervals and the desorption barrier is at 0.57 eV. 
The extra electron in the anti-bonding 2n orbital is seen to 
stretch th C-O bond. The center of mass is moved slightly 
out from the surface in spite of the attraction to the image 
charge. 



tem becomes large when the electronic temperature is 
increased. We also show the probabilities resulting from 
Langevin dynamics with classical and quasiclassical ini- 
tial conditions and it is again seen that the classical ini- 
tial conditions severely underestimates the probabilities. 
In contrast to the harmonic oscillator, the quasiclassical 
approach is in very good approximation for the quan- 
tum initial conditions when calculating desorption prob- 
abilities. This is due to the fact that the quasiclassical 
approach is a good approximation for average quantities 
and the desorption probability in the present case is an 
integral over a continuum of excited states |fc). This will 
be extremely useful since the quasiclassical approxima- 
tion allows us to circumvent phase space sampling. 



IV. AB INITIO POTENTIAL 

As an example illustrating quantum effects in Langevin 
dynamics using ab initio potentials, we consider CO ad- 
sorbed on Cu(lOO). This system has previously been in- 
vestigated in the context of electronic friction and the 
closely connected vibrational linewidth broade ning in- 
duced by electron hole pair excitations, [ll, 113, l35j | All 
parameters in the model Hamiltonian ([1} was obtained 
within Density Functional Theory (DFT) using the code 
gpaw. [36i [37} which is a real-space Density Functional 
Theory (DFT) code that uses the projector augmented 
wave method. m. 111] We used a grid spacing of 0.2 A 
and the calculations were performed in a (2x2) supercell 
sampled by a (4x6) grid of fc-points using the RPBEpoj 
exchange correlation functional. The system was mod- 
elled by a three layer Cu(lOO) slab with the top layer 
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FIG. 6: Density of states projected onto the 2tt orbital of 
CO adsorbed on Cu(lOO) top site. The fuU width at half 
maximum is estimated to be F = 2.0 eV. The Fermi level is at 
E = eV and the resonance is seen to be mostly unoccupied 
in the electronic ground state. 




z-zo [A] 



relaxed and CO adsorbed in a c(2x2) structure (0.5 cov- 
erage at top sites) . For this system the electronic friction 
is dominated by the unoccupied 27r orbitals which we as- 
sume to represent the resonant state |a). 

We have calculated the potential energy surfaces in 
terms of the center of mass and bond length coordi- 
nates which are denoted by z and d respectively. We 
restrict the analysis to these modes since in a first order 
Taylor expansion of ea{xi), the frustrated rotations and 
translations do not couple to the resonant state due to 
symmetry. The desorption energy is determined to be 
Edes ^ 0.57 eV in excellent agreement with the experi- 
mental value 0] ■ The excited state potential energy sur- 
face Vi((i, z) was calculated using a generalization of the 
A-self-consistent field method where the resonant state 
is expanded in a basis of Kohn-Sham orbitals and the 
resulting resonant density is added to the density in each 
iteration step. Thus for each adsorbate position we cal- 
culate the energy resulting from forcing an electron into a 
27r orbital which is then not an eigenstate of the full elec- 
tronic system. The excited state thus has a finite lifetime 
which in the wide band limit can be related to the reso- 
nance width as T = ?i/r.[l^ For details on the method 
and comparison with experiments we refer to. fjlj Since 
electrostatic interactions may arise between an excited 
molecule and its periodic image we have checked that 
the excited state calculations do not change significantly 
when the supercell is changed to (4x4). 

The ground and excited state potential energy surfaces 
are shown in figure [5] The ground state is well approxi- 
mated by a quadratic potential in the internal mode and 
a Morse potential in the center of mass mode. The two 
modes are nearly decoupled and in Tab. |T]we display the 
parameters associated with the two modes at the ground 
state minimum. The resonance width F was obtained 
from the projected density of states shown in figure [51 
At the ground state equilibrium position the width is 
approximately Fq « 2 eV and varying the adsorbate po- 
sition shows that the coordinate dependence is well ap- 
proximated by F = Foe~^/^^ with zr ~ 0.7 A. Since the 
friction tensor is additive in contributing orbitals, we can 



FIG. 7: Diagonal components of the friction tensor as a func- 
tion of COM distance to surface evaluated at T = 6000 K. 
Both components decrease exponentially far from the surface 
but have very different behavior near the minimum position. 



simply multiply the expression by a factor of four to 
account for the degeneracy of the 27r orbital and spin, or 
equivalently, multiply the frictional force by a factor of 
two which for the internal mode reproduces the parame- 
ters used in Sec. IIIII The excitation energy at the ground 
state minimum is Sq = 2.6 eV. The diagonal elements 
of the friction tensor Eq. ^ at the equilibrium position 
and zero temperature can be roughly related to the vi- 
brational lifetimes of the modes: Ti = Mi/rja. In Fig. [7] 
we show the two diagonal components as a function of 
distance to the surface. The two components have the 
same order of magnitude near the equilibrium position 
(z — zq = 0), but the friction in the internal mode (jy^d) is 
seen to decay much faster far from surface than the COM 
friction. Furthermore, the COM friction has a local max- 
imum beyond the equilibrium position and the molecule 
is thus likely to dissipate energy on the path leading to 
desorption which decreases the desorption probability. It 
should be noted that although the frictional force param- 
eters fi have the same order of magnitude, they originate 
from different terms in Eq. (|10p . The center of mass min- 
imum is nearly unaffected by a transition to the excited 
state as seen in figure [5] and the frictional force arises 
only from the COM dependence of the resonance width. 
On the other hand, the resonance width is nearly inde- 
pendent of the internal stretch mode and the internal 
frictional force originate in the large displacement of the 
excited state minimum position. The vibrational life- 
times are in good agreement wi th p revious calculations 
using a different method, [li, S, Eg 

To model a particular surface experiment where a fem- 
tosecond laser pulse induces a surface reaction, one would 
need a detailed model for the time dependent distribu- 
tion of hot electrons resulting from the laser pulse. In 
the present paper we do not aim at a precise quantitative 
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M,/rju{0;0) 


Internal 


0.248 eV 


4.3 eV/A 


2.7 ps 


COM 


0.043 eV 


-3.6 eV/A 


16 ps 



TABLE I: Parameters for the internal vibration and center of 
mass mode for CO adsorbed on a Cu(lOO) top site. 



calculation of reaction rates, but rather wish to examine 
the qualitative impact of including quantum initial states 
in the dynamics. Therefore, we will take a very simple 
model for the hot electrons and assume a thermal pulse 
with a Gaussian temporal shape Te{t) — Tmaxe^'^ /2At- 
with Tjnax = 4000 K and Ai = 0.5 ps. Under the influ- 
ence of this pulse we have performed Langevin dynamics 
with classical quasiclassical and quantized initial condi- 
tions in both the internal and center of mass mode using 
the potentials shown in Fig. [5l The Langevin equation 
is integrated from 2 ps prior to the center of the pulse to 
4 ps after the center of the pulse. Due to the very weak 
coupling between the two modes the initial condition of 
the internal mode has almost no influence on desorption 
probabilities. With fully quantized initial conditions (vi- 
brational ground state) of the COM mode we find a des- 
orption probability of Pquan — 3.7 x 10~^, whereas we 
find Pqc < xlO~^ and Pcias < 10~^ when using quasi- 
classical and classical initial conditions respectively (10® 
trajectories did not result in a single desorption event). 
We note, that when calculating the fluctuating forces Eq. 
dill) , it is most important to take into account the cor- 
relation between the two modes determined by the off- 
diagonal elements of the friction tensor. 

Although a quantization of the internal mode does not 
influence the desorption probability it may have a large 
impact on the distribution of vibrational states of the des- 
orbed molecules. This is illustrated in Fig. [5] where the 
distribution of energy is shown for desorbed molecules 
using the classical, quasiclassical, and quantum initial 
conditions. Due to the low desorption probabilities we 
had to start the molecule with a COM momentum of 
p — SpQ corresponding to 0.19 eV, since otherwise we 
were not able to get good statistics for the energy distri- 
bution of desorbed molecules. However, because of the 
very weak coupling between the two modes, we do not ex- 
pect this to have a large influence on the internal energy 
distribution. The COM energy is not influenced by the 
initial conditions in the internal mode and the difference 
in total energy distributions is solely due to differences in 
the internal mode distributions. It is seen that the qua- 
siclassical initial conditions yields a distribution which 
is similar to the quantized initial conditions, but with 
slightly more weight at high lying energies. The classical 
initial conditions yields a distribution which is inconsis- 
tent with a quantized picture, since from Eq. (jl9l) it 
follows that dP/dE{E = 0) < E^;^ - 8 eV'^. 

To see this in more detail we calculate the probabili- 
ties of the desorbed molecules being in a particular vibra- 
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FIG. 8: The differential probability of desorbed molecule hav- 
ing a given amount of energy as a result of a Gaussian pulse 
of hot electrons with Tmax = 4000 K obtained with classical, 
quasiclassical, and quantized initial conditions. Left is the 
vibrational energy and right is the total energy. 

tional state using the method of Sec. |lll] and Eq. (fTTl) . 

The classical initial conditions lead to po > 1 and pi < 
whereas quasiclassical initial conditions givepi/po — 0.22 
and quantized initial conditions give pi/po — 0.092 which 
is in agreement with Ref. 7. In general, quasiclassical ini- 
tial conditions tend to overestimate pi and underestimate 
Po and p2 as is seen in Fig. |31 In the present case the 
error on pi/po is more than a factor of two. For long in- 
teraction times and high temperatures the quasiclassical 
approximation becomes better and we have repeated the 
above analysis with Tmax — 6000 K, which yields close 
agreement between the vibrational probabilities resulting 
from quasiclassical and quantized initial conditions. 



V. DISCUSSION 

In section Hill it was shown that in order to obtain the 
correct vibrational probabilities for a harmonic oscilla- 
tor, it is crucial to use quantized initial conditions. How- 
ever, quasiclassical initial conditions yield good results 
for the average energy of the harmonic oscillator as well 
as for the desorption probability of the Morse potential. 
Naturally, the quasiclassical approximation is highly at- 
tractive since it only requires a single initial phase space 
point, whereas the correctly quantized initial conditions 
requires a full phase space sampling. In the present work 
we needed a 6 x 6 grid and 10 x 6 grid of initial phase 
space points to represent the relevant part of phase space 
of the harmonic and Morse potentials respectively and 
quantized initial conditions thus required a factor of 36- 
60 more calculations than the quasiclassical approach. 
In general we expect that average quantities are well de- 
scribed by the quasiclassical initial conditions. Similarly, 
high temperatures (compared to the quantum of oscilla- 
tion) and long timescales tend to justify the quasiclassical 
approach. 

With CO on Cu(lOO) as a generic example of a two- 
dimensional problem with ab initio potentials, we found 
that quantization of the internal mode had almost no 
effect on desorption probabilities. However, this is most 
likely due to the weak coupling between the two modes in 
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the present example, but for reactions with very strong 
couphng between modes such as associative desorption 
processes, [i^ quantization of the internal mode is 
likely to be important. Furthermore, if one is interested 
in the final state distribution of vibrational states, it will 
be crucial to take into account the initial zero point mo- 
tion of the adsorbate. For example, the fact that hot 
electron induced associative desorption yields of Hydro- 
gen from Ru(OOOl) are well described by Langevin dy- 
namics except for too low values of desorbate transla- 
tional energies, [2^ H^l may very well be due to initial 
zero point motion. 

It should be mentioned that it is also possible to cal- 
culate the friction tensor directly from density functional 
theory using a basis of Kohn-Sham orbitals. |25l - t27| While 
that method is probably more accurate, the present ap- 
proach based on the reduced density matrix and Newns- 
Anderson like Hamiltonian Eq. ([T|), gives better access 
to the physics involved. For example, in the Newns- 
Anderson framework it is evident that the frictional 
forces on the center of mass mode and the internal mode 
have very different physical origins. On the other hand, 
since the Kohn-Sham approach does not make any as- 
sumption about the physical nature of the friction, it will 
automatically include all contributing states and thus 
give better results when multiple adsorbate states con- 
tribute to the friction. The method applied in the present 
paper only takes into account a single resonance which 
we assume to have a Lorentzian shape, but the excitation 
energy is calculated using ASCF which gives a much bet- 
ter description than the Kohn-Sham eigenvalues. pl| At 
low temperatures, however, the friction is dominated by 
the projected density of states at the Fermi level which is 
unlikely to be well described within the wide band limit 
imposed here. 

We have investigated the importance of including the 
quantized initial state in Langevin dynamics where the 
friction and stochastic force originate from a thermal 
bath of hot electrons. In the title we have referred to 
this as quantum corrected Langevin dynamics, but other 
quantum corrections may also be important. In particu- 
lar, for non-quadratic potentials the time evolution is not 
classical and the Langevin equation should be thought of 
as a semiclassical approximation to the true dynamics. 
In principle, the validity of this approximation should 
always be analyzed in detail for a given potential and 
time of propagation, but very often one can use a quick 
'large n' or similar argument to justify the approxima- 
tion. For example, in the case of CO on Cu(lOO) we ex- 
pect the semiclassical approximation to work well, since 
the Morse potential describing the desorption coordinate 
has 27 bound states within the 0.57 eV potential well, 
which gives an energy spacing much smaller than the av- 
erage adsorbate energy. 

Another quantum effect is that of memory in the fiuc- 
tuating forces. The Markov approximation leading to Eq. 
PTjl completely neglects any correlation between forces 
at different times and essentially only contains thermal 



fluctuations. That the Markov approximation has a clas- 
sical flavor can be seen in the low temperature limit where 
the fluctuating forces vanish. The Langevin equation 
with a harmonic potential then gives rise to a decaying 
average energy: E{t) = i^oe"''*/*^ which is not allowed 
quantum mechanically, since the average energy can not 
become less than the zero point energy. This paradox is 
solved by going beyond the Markov approximation where 
a small fluctuating force exactly cancels the frictional de- 
cay. To get an idea of the range of temperatures where 
the Markov approximation works, we can estimate the 
correlation time by tc — h/ksT . [l9| The timestep used 
in the molecular dynamics in this work was 1 f s which 
corresponds to T = 2900 K and this gives an estimate 
on the lower temperature limit to the Markov approxi- 
mation. Memory effects in non-adiabatic dynamics will 
be explored further in a future paper. 



VI. SUMMARY 

We have analyzed the effect of including zero point mo- 
tion properly in Langevin dynamics with a temperature 
dependent friction tensor. The method which involves 
initial phase space sampling, have been compared to a 
quasiclassical approach where classical initial conditions 
matching the zero point energy is used. For a harmonic 
oscillator, the initial conditions is the only quantum me- 
chanical correction since the quantum dynamics becomes 
classical and we have shown how to obtain vibrational 
probabilities from the classical energy distribution result- 
ing from Langevin dynamics with phase space sampling. 
As expected, the result agrees extremely well with an in- 
herently quantum mechanical master equation approach 
when the initial conditions is included correctly, whereas 
the quasiclassical approach only tends to a reasonable re- 
sult after ^ 1 ps of interaction. We have also compared 
the results of using quantized and quasiclassical initial 
conditions in a Morse potential and found only a little 
effect on the probability for escaping the potential well. 
The reason for this can be attributed to the fact that the 
escape probabilities involves a sum over the continuous 
set of free states and the quasiclassical approach yields a 
good description of average quantities. 

With CO on Cu(lOO) as a generic example, we have 
demonstrated the effect for an adsorbate system with ab 
initio potentials and electronic friction. The tensor struc- 
ture of the friction introduces correlation in the fluctuat- 
ing forces and the nonlinear interaction gives rise to posi- 
tion dependence in the friction. For a model pulse of hot 
electrons we showed that, compared to the quasiclassi- 
cal approach, quantized initial conditions both increases 
the desorption probability and changes the distribution 
of vibrational states. 
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Appendix A: Path integral representation of the 
harmonic oscillator density matrix 

The derivation of Langcvin dynamics is most easily 
done with a path integral representation of the reduced 
time dependent density matrix. To see how it works we 
consider again the harmonic oscillator with a simple de- 
gree of freedom and start with the expression ([S]). The 
two propagators can be written as path integrals result- 
ing in 



p{x,y;t) ^ J dxodyQ{xa\pa\yo) (Al) 

with the action 

So[x{t')] = J\t'(^mx^t') - ^muj'x'it')), (A2) 

and x{0) — xq, j/(0) ~ yo. Introducing the average path 
u{t) = x(t)/2 -\- y{t)/2 and the fluctuation v{t) = x{t) — 
y{t) we can do a partial integration on the kinetic term 
and write the sum of actions 

S'o[a;(t')] — SQ[y{t')] — j dt' (^miiv — muj'^uvj 

dt' (^mil + muj^uj v. 



— muv 

"'0 

Thus the density matrix becomes 

pix,y;t) = / dxodyoixolpolyo) 



(A3) 



It is now straightforward to perform the path integral in 
v{t') which gives a delta functional on the classical path 
■u{t') — —uPu{t') for the average coordinate. If we are 
only interested in the probabilities of finding the particle 
at a given position we just need the diagonal elements of 
the density matrix where the end points satisfy u{t) — 
x{t) = y{t) and v{t) = 0. In terms of these coordinates 
the diagonal part of the density matrix becomes 



p{u;t) (X J duodvo{uo + vo/2\po\uQ-vo/2) 

(X J duoV{uo,po{uo,u{t))), (A4) 



where 



r{x,p) = ^J dy{x + y/2\po\x ~ yl2)e 



■ipy/h 



(A5) 



is the Wigner distribution of an initial state described by 
the density matrix pQ and the path integral delta function 
has been collapsed by noting that for a given uq there is 
a unique initial momentum po given by 



Po = muo 



muj 



sin ujt 



{u(t) — Uo COSLOt) 



(A6) 



that connects the initial position classically with u{t). 
The easiest way to determine the normalization is to re- 
quire that J dup{u;t) ~ 1, and the expression is then 
seen to be identical to Eq. ([6]). 
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